1 Figures
1.1 Map
knitr::include_graphics('figs/statmap.jpg')
Figure 1.1: Locations of stations in the Puget Sound where pteropod and environmental sampling occurred. Samples were collected in April, July, and September from 2014 to 2016.
1.2 Site pteropod distributions
colmos <- diverging_hcl(3, 'Blue-Red')
# make cohortyr an ordered factor
toplo <- biodat %>%
mutate(cohortyr = factor(cohortyr, ordered = T)) %>%
select(mo, station, yr,typ3, len) %>%
gather('var', 'val', typ3, len) %>%
mutate(
var = factor(var, levels = c('len', 'typ3'), labels = c('length (um)', '% type III dissolution'))
)
p <- ggplot(toplo, aes(x = mo, y = val, fill = )) +
geom_boxplot(aes(fill = mo), outlier.shape = NA, alpha = 0.4) +
geom_label(aes(fill = mo, label = station), size = 4, alpha = 0.6, position = position_jitter(width = 0.3),
label.padding = unit(0.1, "lines")) +
facet_wrap(~var, ncol = 2, scales = 'free_y', strip.position = 'left') +
scale_fill_manual(values = colmos) +
# scale_colour_manual(values = colmos) +
theme_bw(base_family = 'serif', base_size = 14) +
theme(
strip.placement = 'outside',
axis.title = element_blank(),
strip.background = element_blank(),
legend.position = 'none'
)
jpeg('figs/biodst.jpg', height = 5, width = 8, units = 'in', res = 300, family = 'serif')
print(p)
dev.off()knitr::include_graphics('figs/biodst.jpg')
Figure 1.2: Pteropod distributions by length (left) and dissolution (right) across sites and months. Site numbers are shown in each point. Box plots represent the median values with the upper and lower limits of the boxes defined as the 25th and 75th percentile of the distributions. The whiskers are 1.5 times the interquartile range.
1.3 Cluster analyses
depth_dat <- read_excel('raw/WOAC PS Cruise Plan.xlsx', sheet = 'Sheet2') %>%
select(station = Station, depth = `Depth (m)`) %>%
mutate(
station = gsub('^P', '', station),
station = as.numeric(station)
)
toplo <- read_excel('raw/WOAC_data_5-1-2018_for_Nina.xlsx', sheet = 'ALL_DATA', na = c('', '-999')) %>%
select(Date_collected, STATION_NO, LATITUDE_DEC, LONGITUDE_DEC, NISKIN_NO, CTDTMP_DEG_C_ITS90,
CTDSAL_PSS78, CTDOXY_UMOL_KG_ADJ, `Omega Ar`) %>%
rename(
date = Date_collected,
station = STATION_NO,
lat = LATITUDE_DEC,
lon = LONGITUDE_DEC,
niskin = NISKIN_NO,
temp = CTDTMP_DEG_C_ITS90,
sal = CTDSAL_PSS78,
ara = `Omega Ar`,
oxy = CTDOXY_UMOL_KG_ADJ
) %>%
left_join(depth_dat, by = 'station') %>%
gather('var', 'val', temp:depth) %>%
group_by(station, var) %>%
summarize(
valmn = mean(val, na.rm = T)
) %>%
ungroup %>%
nest(data = everything()) %>%
mutate(
disval = purrr::map(data, function(x){
x %>%
select(station, var, valmn) %>%
spread(var, valmn) %>%
data.frame %>%
column_to_rownames('station') %>%
decostand(method = 'standardize') %>%
vegdist(method = 'euclidean')
}),
clsval = purrr::map(disval, function(x){
x %>%
hclust(method = 'average')
}),
cutval = purrr::map(clsval, function(x){
# get cut groups
cutree(x, k = 3)
}),
denplo = purrr::pmap(list(cutval, clsval), function(cutval, clsval){
# get order
clstord <- order.hclust(clsval) %>%
cutval[.] %>%
unique
# get colors, correct by order
cols <- colmst %>%
.[clstord]
pdend <- clsval %>%
as.dendrogram %>%
set("branches_k_color", k = 3, value = cols) %>%
set("labels_colors", k = 3, value = cols) %>%
set("labels_cex", 0.8)
p1 <- as.ggdend(pdend) %>%
ggplot(horiz = TRUE, offset_labels = -0.1)
p1
}),
displo = purrr::pmap(list(disval, cutval, clsval), function(disval, cutval, clsval){
# prep distance data to plot
# long format of dist matrix
toplo <- disval %>%
as.matrix %>%
as.data.frame %>%
rownames_to_column('station') %>%
gather('station2', 'dist', -station) %>%
arrange(dist) %>%
mutate(
dist = ifelse(station == station2, NA, dist)
)
# get site order levels based on clustering
sitfc <- clsval$labels[clsval$order]
toplo <- toplo %>%
mutate(
station = factor(station, levels = sitfc),
station2 = factor(station2, levels = sitfc)
)
# plot
p <- ggplot(toplo) +
geom_tile(aes(x = station, y = station2, fill = dist), colour = 'black') +
scale_x_discrete('', expand = c(0, 0)) +
scale_y_discrete('', expand = c(0, 0)) +
scale_fill_gradient2('Dissimilarity between stations\nby water chemistry', low = 'lightblue', mid = 'white', high = 'tomato1', midpoint = 2.5, limits = c(0.5, 4.5)) +
guides(fill = guide_colourbar(barheight = 0.5, barwidth = 10, label.theme = element_text(size = 11, angle = 0))) +
pbase
# index values of cluster cuts
brks <- cutval %>%
.[clsval$order] %>%
duplicated %>%
`!` %>%
which %>%
`-` (0.5) %>%
.[-1]
p <- p +
geom_vline(xintercept = brks, size = 1.5) +
geom_hline(yintercept = brks, size = 1.5) +
theme_bw(base_family = 'serif') +
theme(
legend.position = 'top',
legend.direction = 'horizontal'
)
return(p)
}),
disleg = purrr::map(displo, function(x) g_legend(x)),
displo = purrr::map(displo, function(x){
p <- x + theme(legend.position = 'none')
return(p)
}),
mapplo = purrr::pmap(list(clsval, cutval), function(clsval, cutval){
mapplo <- locs %>%
mutate(cutval = rev(cutval))
# plot the basemap
p <- ggmap(bsmap) +
geom_point(data = mapplo, aes(x = lon, y = lat, fill = factor(cutval)), pch = 21, size = 8, alpha = 0.7) +
geom_text(data = mapplo, aes(x = lon, y = lat, label = station), colour = 'white', size = 3) +
scale_fill_manual(values = colmst) +
theme_bw(base_family = 'serif') +
theme(
axis.title = element_blank(),
legend.position = 'none',
axis.text.y = element_text(size = 8),
axis.text.x = element_text(size = 8, angle = 45, hjust = 1)
)
return(p)
})
)
jpeg('figs/clsts.jpg', height = 3.5, width = 7, units = 'in', res = 300, family = 'serif')
wrap_elements(toplo$disleg[[1]]) + (toplo$denplo[[1]]+ toplo$displo[[1]] + toplo$mapplo[[1]])+ plot_layout(ncol = 1, heights = c(0.05, 1))
dev.off()knitr::include_graphics('figs/clsts.jpg')
Figure 1.3: Clustering results of stations based on annual averages of salinity, water temperature, dissolved oxygen, and aragonite saturation state. Station-depth is also included. Averages are based on all environmental data collected across the sample years from 2014 to 2016 in the same month. Results are shown as dendrograms for site clustering (left), dissimilarity matrices showing mean Euclidean distances between observations at pairs of sites (middle), and spatial arrangements of the defined clusters (right). Cluster groups were set at three based on approximate dendrogram separation between sites to explain dominant patterns among environmental variables.
1.4 PCA
depth_dat <- read_excel('raw/WOAC PS Cruise Plan.xlsx', sheet = 'Sheet2') %>%
select(station = Station, depth = `Depth (m)`) %>%
mutate(
station = gsub('^P', '', station),
station = as.numeric(station)
)
# pteropod birthday
strt <- '2008-06-01' %>%
as.Date
# make cohortyr an ordered factor
biodat <- biodat %>%
mutate(cohortyr = factor(cohortyr, ordered = T))
chmdatsum <- chmdatsum %>%
mutate(cohortyr = factor(cohortyr, ordered = T))
# combine data for pca
biosub <- biodat %>%
select(cohortyr, mo, station, typ3)
chmsub <- chmdatsum %>%
filter(var %in% c('ara', 'temp', 'oxy', 'sal')) %>%
select(-date, -yr, -lon, -lat, -max, -std, -rng, -dlt) %>%
gather('valtyp', 'val', ave, min) %>%
filter(var %in% 'ara' & valtyp %in% 'min' | !var %in% 'ara' & valtyp %in% 'ave') %>%
select(-valtyp) %>%
spread(var, val)
tomod <- chmsub %>%
left_join(biosub, by = c('cohortyr', 'mo', 'station')) %>%
left_join(clsts, by = 'station') %>%
left_join(depth_dat, by = 'station') %>%
unite('stat_mo', station, mo, sep = ', ', remove = F) %>%
unite(stat_moyr, stat_mo, cohortyr, sep = ' ', remove = F) %>%
filter(!is.na(typ3)) %>%
as.data.frame(stringsAsFactors = F) %>%
column_to_rownames('stat_moyr')
# pc mod
mod <- prcomp(tomod[ , c('ara', 'oxy', 'temp', 'sal', 'depth')], scale. = T, center = T)
# no labels
p1 <- ggord(mod, grp_in = as.character(tomod$clst), vec_ext = 4, size = tomod$typ3, coord_fix = F, labcol = 'blue', alpha = 0.85) +
scale_size(range = c(2, 8)) +
scale_colour_manual(values = rev(colmst)) +
scale_fill_manual(values = rev(colmst)) +
guides(size = guide_legend(title = '% type III dissolution')) +
theme(legend.position = 'top')
pleg <- g_legend(p1)
p1 <- p1 + theme(legend.position = 'none')
p2 <- ggord(mod, axes = c('1', '3'), grp_in = as.character(tomod$clst), vec_ext = 4, size = tomod$typ3, coord_fix = F, labcol = 'blue', alpha = 0.85) +
scale_size(range = c(2, 8)) +
scale_colour_manual(values = rev(colmst)) +
scale_fill_manual(values = rev(colmst)) +
guides(size = guide_legend(title = '% type III dissolution')) +
theme(legend.position = 'none')
# with labels
# p2 <- ggord(mod, obslab = T, vec_ext = 4, size = 1.8, coord_fix = F, labcol = 'blue')
jpeg('figs/pcastat.jpg', family = 'serif', height = 4.5, width = 8.5, res = 300, units = 'in')
grid.arrange(
pleg,
arrangeGrob(p1, p2, ncol = 2),
ncol = 1, heights = c(0.1, 1)
)
dev.off()knitr::include_graphics('figs/pcastat.jpg')
Figure 1.4: Results of principal components analysis for environmental variables collected at each site for each sample date. Environmental variables included temperature, salinity, dissolved oxygen, and minimum aragonite saturation state. Station depth is also included. The left plot shows site groupings based on dominant clusters shown in Figure 1.3, with site points sized by measured type III dissolution for pteropods collected at the same location and date. The right plot shows the same as the first plot but for the first and third principal component axes.
1.5 CTD profiles
# cluster
clsts <- tibble(
station = c(402, 38, 28, 22, 12, 8, 4),
clst = c(3, 2, 2, 1, 3, 2, 3)
)
# date levels, labels
dtlev <- as.character(sort(unique(chmdatraw$date)))
dtlab <- as.character(format(unique(chmdatraw$date), '%Y-%m'))
# chemistry raw data
toplo <- chmdatraw %>%
filter(var %in% 'ara') %>%
left_join(clsts, by = 'station') %>%
filter(!is.na(val)) %>%
group_by(date, mo, station, clst, depth, var) %>%
summarise(val = mean(val)) %>%
ungroup %>%
mutate(
date = factor(date, levels = dtlev, labels = dtlab),
alph = case_when(
val > 1.2 ~ 0.3,
val > 1 & val <= 1.2 ~ 0.3,
val <= 1 ~ 1
)
) %>%
na.omit
p <- ggplot(toplo, aes(fill = factor(clst), x = val, y = depth)) +
geom_point(pch = 21, colour = 'black', size = 3, alpha = toplo$alph) +
# geom_smooth(method="nls",
# formula = y ~ SSasymp(x, yf, y0, log_alpha),
# se=FALSE) +
scale_fill_manual(values = rev(colmst)) +
geom_vline(linetype = 'dashed', aes(xintercept = 1), size = 1) +
geom_vline(linetype = 'dotted', aes(xintercept = 1.2), size = 1) +
scale_y_reverse() +
scale_x_continuous(breaks = seq(0, 5, by = 0.5)) +
ylab('Depth (m)') +
xlab(expression(Omega['ar,min'])) +
facet_grid(clst~mo) +
theme_bw(base_family = 'serif', base_size = 14) +
theme(
strip.background = element_blank(),
legend.title = element_blank(),
legend.position = 'none',
# panel.grid = element_blank(),
panel.grid.minor = element_blank()#,
# axis.text.x = element_text(size = 8)
)
jpeg('figs/ctdplo.jpg', height = 6.5, width = 8, units = 'in', res = 300, family = 'serif')
p
dev.off()knitr::include_graphics('figs/ctdplo.jpg')
Figure 1.5: Depth profiles of minimum aragonite saturation state across seasons and stations. Stations are grouped by clusters in Figure @fig(fig:clstmap). Vertical lines indicate minimum saturation state of 1 (dashed) and 1.2 (dotted). Point transparencies also correspond to observations above or below minimum saturation state of 1.
1.6 Observed time series
strdat <- chmdatsum %>%
filter(var %in% 'ara') %>%
full_join(biodat, by = c('station', 'date', 'yr', 'cohortyr', 'mo')) %>%
select(-lon, -lat, -var, -abu, -len, -yr) %>%
gather('ara', 'chmval', ave:dlt) %>%
filter(ara %in% 'min') %>%
mutate(
thrsh = 1.2,
wts = case_when(
mo %in% c('Jul', 'Sep') ~ 1,
mo == 'Apr' ~ 3.5
),
aracat = ifelse(chmval < thrsh, 1, 0),
aradff = thrsh - chmval
) %>%
unite('costa', cohortyr, station, remove = F) %>%
group_by(costa) %>%
mutate(
strsdis = cumsum(aracat),
strscnt = cumsum(aradff),
strscntwt = cumsum(aradff * wts)
) %>%
ungroup %>%
mutate(
strscntwt = scales::rescale(strscntwt, to = c(0, 1))
) %>%
left_join(clsts, by = 'station') %>%
unite('station', station, clst, sep = ' (', remove = F) %>%
mutate(station = paste0(station, ')'))
# bio
toplo1 <- biodat %>%
select(date, yr, cohortyr, mo, station, typ1, typ2, typ3) %>%
gather('diss', 'val', typ1:typ3) %>%
filter(diss %in% 'typ3') %>%
left_join(clsts, by = 'station') %>%
unite('station', station, clst, sep = ' (', remove = F) %>%
mutate(station = paste0(station, ')'))
p1 <- ggplot(toplo1, aes(x = factor(date), y = val)) +
geom_line(aes(group = cohortyr)) +
geom_point(aes(fill = mo), alpha = 0.7, size = 3, pch = 21) +
facet_grid(reorder(station, clst)~.) +
theme_bw(base_family = 'serif', base_size = 10) +
theme(
strip.background = element_blank(),
axis.title.x = element_blank(),
legend.title = element_blank(),
axis.text.x = element_text(size = 8, angle = 45, hjust = 1),
# panel.grid = element_blank(),
legend.position = 'none'
) +
ylab('% type III dissolution') +
scale_x_discrete(labels = format(sort(unique(toplo1$date)), '%Y-%m')) +
scale_fill_manual(values = colmos)
# chem
toplo2 <- chmdatsum %>%
filter(var %in% 'ara') %>%
left_join(clsts, by = 'station') %>%
unite('station', station, clst, sep = ' (',remove = F) %>%
mutate(station = paste0(station, ')'))
p2 <- ggplot(toplo2, aes(x = factor(date), y = min)) +
geom_segment(data = strdat, aes(y = 1.2, yend = chmval, xend = factor(date), colour = 'Cumulative stress difference'), linetype = 'dotted', size = 0.9) +
geom_line(aes(group = cohortyr)) +
geom_point(aes(fill = mo), alpha = 0.7, size = 3, pch = 21) +
facet_grid(reorder(station, clst)~.) +
theme_bw(base_family = 'serif', base_size = 10) +
theme(
strip.background = element_blank(),
axis.title.x = element_blank(),
# legend.title = element_blank(),
axis.text.x = element_text(size = 8, angle = 45, hjust = 1),
# panel.grid = element_blank(),
legend.position = 'top'
) +
guides(guide_legend) +
scale_color_discrete(name="") +
scale_fill_manual('Cohort months', values = colmos) +
ylab(expression(Omega['ar,min'])) +
scale_x_discrete(labels = format(sort(unique(toplo2$date)), '%Y-%m')) +
geom_text(data = strdat, aes(y = 1, label = round(strscntwt, 1)), colour = 'tomato1', vjust = 1.2, size = 3, hjust = -0.2) +
geom_hline(aes(yintercept = 1.2), colour = 'tomato1', size = 1, alpha = 0.7)
pleg <- g_legend(p2)
p2 <- p2 + theme(legend.position = 'none')
jpeg('figs/obsdat.jpg', height = 6, width = 7, units = 'in', res = 300, family = 'serif')
wrap_elements(pleg) +
(p1 + p2 + plot_layout(ncol = 2)) +
plot_layout(ncol = 1, heights = c(0.05, 1))
dev.off()knitr::include_graphics('figs/obsdat.jpg')
Figure 1.6: Observed time series for each station (rows, group number from Figure 1.3 in parentheses) showing % type III dissolution of pteropods (left) and observed minimum aragonite saturation state (right). Points at each station are connected by cohort years. The right plot shows the selected aragonite threshold (\(\Omega_{crit} = 1\)) as a horizontal line with the difference between the threshold and minimum observed value shown as a dotted line.
1.7 Dissolution vs aragonite
toplo <- chmdatsum %>%
filter(var %in% 'ara') %>%
full_join(biodat, by = c('station', 'date', 'yr', 'cohortyr', 'mo')) %>%
select(-lon, -lat, -var, -abu, -len) %>%
gather('dissvar', 'dissval', typ1:typ3) %>%
gather('chemvar', 'chemval', ave:dlt) %>%
filter(dissvar %in% 'typ3') %>%
filter(chemvar %in% 'min')
p1 <- ggplot(toplo, aes(x = chemval, y = dissval)) +
# geom_line(aes(group = station), colour = 'grey') +
stat_smooth(method = 'lm', se = T, colour = 'black') +
geom_point(aes(colour = factor(cohortyr), group = station), pch = 15, size = 0, alpha = 0) +
geom_label(aes(fill = factor(cohortyr), group = station, label = station), colour = 'black', size = 3, alpha = 0.6, show.legend = F) +
facet_grid(~ mo) + #, scales = 'free') +
theme_bw(base_family = 'serif', base_size = 12) +
theme(
strip.background = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
axis.title = element_blank()
) +
scale_fill_manual(values = colmos) +
scale_colour_manual(values = colmos) +
guides(colour = guide_legend(override.aes = list(size = 4, alpha = 0.6)))
p2 <- ggplot(toplo, aes(x = chemval, y = dissval)) +
# geom_line(aes(group = station), colour = 'grey') +
stat_smooth(method = 'lm', se = T, colour = 'black') +
geom_point(aes(colour = factor(mo), group = station), pch = 15, size = 0, alpha = 0) +
geom_label(aes(fill = factor(mo), group = station, label = station), colour = 'black', size = 3, alpha = 0.6, show.legend = F) +
facet_grid(~ cohortyr) + #, scales = 'free') +
theme_bw(base_family = 'serif', base_size = 12) +
theme(
strip.background = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
axis.title = element_blank()
) +
scale_fill_manual(values = colmos) +
scale_colour_manual(values = colmos) +
guides(colour = guide_legend(override.aes = list(size = 4, alpha = 0.6)))
jpeg('figs/disvara.jpg', height = 6, width = 7, units = 'in', res = 300, family = 'serif')
wrap_elements(textGrob(' % type III dissolution', rot = 90)) + (p1 + p2 + wrap_elements(textGrob(expression(Omega['ar, min']))) + plot_layout(ncol = 1, heights = c(1, 1, 0.1))) + plot_layout(ncol = 2, widths = c(0.1, 1))
dev.off()knitr::include_graphics('figs/disvara.jpg')
Figure 1.7: Percent type III dissolution measured in pteropods versus minimum observed aragonite saturation state for each station. The top row shows stations grouped by month across cohort years and the bottom row shows stations grouped by cohort years across months. Linear regression lines with 95% confidence intervals are shown in each panel.
Linear model by years.
modyrs <- lm(dissval ~ chemval * factor(cohortyr), data = toplo)#$toplo[!toplo$cohortyr %in% '2016', ])
summary(modyrs)##
## Call:
## lm(formula = dissval ~ chemval * factor(cohortyr), data = toplo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50.725 -12.335 2.017 13.162 34.678
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.545 7.468 6.232 1.54e-07 ***
## chemval -32.222 9.991 -3.225 0.00238 **
## factor(cohortyr).L 3.104 13.536 0.229 0.81968
## factor(cohortyr).Q -16.192 12.306 -1.316 0.19509
## chemval:factor(cohortyr).L -1.346 18.223 -0.074 0.94147
## chemval:factor(cohortyr).Q 8.754 16.336 0.536 0.59474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.48 on 44 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.2703, Adjusted R-squared: 0.1874
## F-statistic: 3.26 on 5 and 44 DF, p-value: 0.01369Linear model by months.
modmos <- lm(dissval ~ chemval * factor(mo), data = toplo)#$toplo[!toplo$cohortyr %in% '2016', ])
summary(modmos)##
## Call:
## lm(formula = dissval ~ chemval * factor(mo), data = toplo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.357 -12.690 2.726 12.136 38.841
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 54.99 8.62 6.380 9.31e-08 ***
## chemval -40.56 11.80 -3.436 0.0013 **
## factor(mo).L 26.50 14.78 1.793 0.0799 .
## factor(mo).Q 15.52 15.08 1.029 0.3090
## chemval:factor(mo).L -19.25 18.38 -1.047 0.3007
## chemval:factor(mo).Q -11.23 22.32 -0.503 0.6175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.88 on 44 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.3144, Adjusted R-squared: 0.2364
## F-statistic: 4.035 on 5 and 44 DF, p-value: 0.00423Combined model, average aragonite.
#data
tomod <- chmdatsum %>%
filter(var %in% 'ara') %>%
full_join(biodat, by = c('station', 'date', 'yr', 'cohortyr', 'mo')) %>%
select(-lon, -lat, -var, -abu, -len) %>%
gather('dissvar', 'dissval', typ1:typ3) %>%
gather('chemvar', 'chemval', ave:dlt) %>%
filter(dissvar %in% 'typ3') %>%
filter(chemvar %in% 'ave')
# model
modall <- lm(dissval ~ chemval, data = tomod)
summary(modall)##
## Call:
## lm(formula = dissval ~ chemval, data = tomod)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.503 -16.234 -2.162 16.564 56.050
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.22 13.46 3.658 0.000631 ***
## chemval -23.87 12.73 -1.875 0.066859 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.07 on 48 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.06826, Adjusted R-squared: 0.04884
## F-statistic: 3.516 on 1 and 48 DF, p-value: 0.06686# plot
plot(dissval ~ chemval, data = tomod)
abline(reg = modall)
# predicted
ave <- mean(tomod$chemval, na.rm = T)
toprd <- data.frame(
chemval = c(ave + 0.4, ave),
lab = c('pre (1750)', 'current')
) %>%
mutate(
type3 = predict(modall, newdata = ., se.fit = T)$fit,
sefit = predict(modall, newdata = ., se.fit = T)$se.fit
) %>%
rename(ara = chemval)
toprd## ara lab type3 sefit
## 1 1.43739 pre (1750) 14.91337 5.969630
## 2 1.03739 current 24.45971 2.980925Combined model, minimum aragonite.
#data
tomod <- chmdatsum %>%
filter(var %in% 'ara') %>%
full_join(biodat, by = c('station', 'date', 'yr', 'cohortyr', 'mo')) %>%
select(-lon, -lat, -var, -abu, -len) %>%
gather('dissvar', 'dissval', typ1:typ3) %>%
gather('chemvar', 'chemval', ave:dlt) %>%
filter(dissvar %in% 'typ3') %>%
filter(chemvar %in% 'min')
# model
modall <- lm(dissval ~ chemval, data = tomod)
# plot
plot(dissval ~ chemval, data = tomod)
abline(reg = modall)
# predicted
ave <- mean(tomod$chemval, na.rm = T)
toprd <- data.frame(
chemval = c(ave + 0.5, ave),
lab = c('pre (1750)', 'current')
) %>%
mutate(
type3 = predict(modall, newdata = ., se.fit = T)$fit,
sefit = predict(modall, newdata = ., se.fit = T)$se.fit
) %>%
rename(ara = chemval)
toprd## ara lab type3 sefit
## 1 1.1956505 pre (1750) 8.352551 5.712529
## 2 0.6956505 current 24.580258 2.7927911.8 Oxygen and aragonite
# input
oxydat <- chmdatsum %>%
filter(var %in% c('ara', 'oxy', 'araund')) %>%
full_join(biodat, by = c('station', 'date', 'yr', 'cohortyr', 'mo')) %>%
select(-lon, -lat, -abu, -len) %>%
gather('dissvar', 'dissval', typ1:typ3) %>%
gather('chemvar', 'chemval', ave:dlt) %>%
filter(dissvar %in% 'typ3')
plo1dat <- oxydat %>%
filter(var %in% 'ara' & chemvar %in% 'min' | var %in% 'oxy' & chemvar %in% 'max') %>%
select(-chemvar) %>%
spread(var, chemval) %>%
filter(mo %in% c('Jul'))
plo1mod <- lm(dissval ~ ara * oxy, data = plo1dat)
# colors
cols <- RColorBrewer::brewer.pal(9, 'RdBu')
# get data to plot
pldat <- get_pldat(plo1mod, 'oxy', fct = 0.75, pos = 'left')
pl1a <- pldat[[1]]
pl1b <- pldat[[2]]
# names for aes_string
nms <- names(pl1a)
p1 <- ggplot() +
geom_ribbon(data = pl1a, aes_string(x = nms[3], ymin = 'lo', ymax = 'hi', group = nms[4]), alpha = 0.5, fill = 'grey') +
geom_line(data = pl1a, aes_string(x = nms[3], y = nms[1], group = nms[4], colour = nms[4]), size = 1) +
geom_text(data = pl1b, aes(x= x, y = y, label = lab), hjust = 0) +
theme_bw() +
scale_x_continuous(expression(Omega['ar, min'])) +
scale_y_continuous('% type III dissolution') +
coord_cartesian(ylim = c(-10, 90), xlim = c(0.19, 0.75)) +
scale_colour_gradientn(expression(paste(O[2], ' (', mu, 'mol ', kg^{-1}, ')')), colours = cols)
jpeg('figs/oxyint.jpg', family = 'serif', height = 4, width = 5.25, res = 300, units = 'in')
p1
dev.off()knitr::include_graphics('figs/oxyint.jpg')
Figure 1.8: Relationships of dissolved oxygen and aragonite saturation state with dissolution measures of pteropods. The plot shows results for a linear model of dissolution against minimum observed aragonite state and maximum oxygen that includes both main and interaction effects (July only, F = 6.1, df = 3,16, R2 = 0.53, p < 0.05). The color range depicts the minimum and maximum observed values for maximum dissolved oxygen across all stations.
summary(plo1mod)##
## Call:
## lm(formula = dissval ~ ara * oxy, data = plo1dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.865 -8.984 -0.347 6.584 24.740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.1257 61.9108 3.168 0.00597 **
## ara -228.2824 75.2017 -3.036 0.00787 **
## oxy -0.4474 0.1769 -2.529 0.02233 *
## ara:oxy 0.5993 0.2275 2.635 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.73 on 16 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.53, Adjusted R-squared: 0.4419
## F-statistic: 6.015 on 3 and 16 DF, p-value: 0.0060561.9 Growth cohorts
toplo <- biodat %>%
mutate(
cdate = ymd(paste(2015, month(date), day(date), sep = '-')),
cohortseasgrp = case_when(
yr == 2014 & cohortseas == 'winter' ~ 'winter2014',
yr == 2015 & cohortseas == 'winter' & mo == 'Apr' ~ 'winter2014',
yr == 2015 & cohortseas == 'summer' ~ 'summer2015',
yr == 2015 & cohortseas == 'winter' & mo == 'Sep' ~ 'winter2015',
yr == 2016 & cohortseas == 'winter' & mo == 'Apr' ~ 'winter2015',
yr == 2016 & cohortseas == 'summer' ~ 'summer2016'
),
cohortseas = factor(cohortseas, levels = c('summer', 'winter'), labels = c('summer', 'fall')),
date = floor_date(date, unit = 'month')
) %>%
filter(!is.na(cohortseasgrp)) %>%
rename(Cohort = cohortseas)
p <- ggplot(toplo, aes(x = date, y = len)) +
geom_point(aes(fill = Cohort), pch = 21, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', aes(group = cohortseasgrp, colour = Cohort, fill = Cohort), show.legend = F) +
geom_text_repel(aes(label = station), size = 2.5) +
theme_bw(base_family = 'serif') +
scale_x_date(date_breaks = '2 months', date_labels = '%b %Y') +
ylab(expression(paste("Length (", mu, "m)"))) +
theme(
axis.title.x = element_blank(),
legend.position = 'top'
)
jpeg('figs/groplo.jpg', family = 'serif', height = 3, width = 7, res = 300, units = 'in')
p
dev.off()knitr::include_graphics('figs/groplo.jpg')
Figure 1.9: Growth cohorts represented as summer and fall individuals with measured lengths against time. Growth trends within each cohort are shown with linear regression lines. Individuals from April observations were assigned to seasonal growth cohorts using a cutoff length of 500 \(\mu\)m between juveniles of the summer cohort and adults of the prior fall cohort. All individuals in July were assumed to be adults of the summer cohort (i.e., no surviving adults from the previous fall cohort). Similarly, all individuals in September were considered juveniles of the fall cohort (i.e., no surviving adults from the previous summer cohort).
1.10 Cumulative stress plot
# rename pteropod response measure to generic
toplo <- strdat %>%
rename(
rsp = 'typ3'
)
p1 <- ggplot(toplo, aes(x = mo, y = strscntwt, group = costa)) +
geom_line(linetype = 'dotted') +
geom_point(aes(colour = mo, group = station), pch = 15, size = 0, alpha = 0) +
geom_point(aes(group = station, size = rsp), pch = 15, alpha = 0, colour = 'black', position = position_jitter(width = 0.1, height = 0.1)) +
geom_label(aes(fill = mo, group = costa, label = station, size = rsp), colour = 'black', position = position_jitter(width = 0.1, height = 0.1), alpha = 0.6, show.legend = F) +
ylab('Cumulative stress magnitude') +
facet_wrap(~cohortyr) +
theme_bw(base_family = 'serif', base_size = 12) +
theme(
axis.title.x = element_blank(),
legend.position = 'top',
strip.background = element_blank()
) +
scale_size('% type III dissolution') +
guides(fill = guide_legend(title = element_blank(), override.aes = list(size = 4))) +
scale_colour_manual(values = colmos) +
scale_fill_manual(values = colmos) +
guides(
colour = guide_legend(title = '', override.aes = list(size = 4, alpha = 0.6)),
size = guide_legend(override.aes = list(size = seq(2, 8, length = 4), alpha = 0.6))
)
p2 <- ggplot(toplo, aes(x = strscntwt, y = rsp)) +
geom_line(aes(group = costa), linetype = 'dotted') +
geom_point(aes(colour = mo, group = station), pch = 15, size = 0, alpha = 0) +
geom_label(aes(fill = mo, group = costa, label = station), size = 3, colour = 'black', alpha = 0.6, show.legend = F) +
geom_smooth(method = 'lm', colour = 'black') +
xlab('Cumulative stress magnitude') +
ylab('% type III dissolution') +
facet_wrap(~cohortyr) +
theme_bw(base_family = 'serif', base_size = 12) +
theme(
legend.position = 'top',
legend.title = element_blank(),
strip.background = element_blank()
) +
scale_colour_manual(values = colmos) +
scale_fill_manual(values = colmos) +
guides(
colour = guide_legend(title = '', override.aes = list(size = 4, alpha = 0.6))
)
jpeg('figs/cumstr.jpg', family = 'serif', res = 300, units = 'in', height = 7, width = 7)
p1 + p2 + plot_layout(ncol = 1, heights = c(1, 1))
dev.off()knitr::include_graphics('figs/cumstr.jpg')
Figure 1.10: Relationships between percent type III dissolution and cumulative stress magnitude within cohort years. The top plot shows the progression of estimated cumulative stress from July to April throughout a cohort year for each station, with points sized by percent dissolution. The bottom plot shows the estimated linear relationship between percent dissolution and cumulative stress. The cumulative stress estimates within a year represent the frequency and magnitude of estimated exposure time of pteropods in a cohort when conditions were under-saturated below \(\Omega_{crit} = 1\).
strmods <- strdat %>%
group_by(cohortyr) %>%
nest %>%
mutate(
mod = purrr::map(data, function(x) summary(lm(typ3 ~ strscntwt, x)))
)
strmods## # A tibble: 3 x 3
## # Groups: cohortyr [3]
## cohortyr data mod
## <ord> <list> <list>
## 1 2014 <tibble [21 x 19]> <smmry.lm>
## 2 2015 <tibble [21 x 19]> <smmry.lm>
## 3 2016 <tibble [14 x 19]> <smmry.lm>strmods$mod## [[1]]
##
## Call:
## lm(formula = typ3 ~ strscntwt, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.212 -12.562 -2.033 12.453 23.418
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.506 6.652 0.978 0.344
## strscntwt 49.742 21.863 2.275 0.038 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.12 on 15 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.2566, Adjusted R-squared: 0.207
## F-statistic: 5.176 on 1 and 15 DF, p-value: 0.038
##
##
## [[2]]
##
## Call:
## lm(formula = typ3 ~ strscntwt, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.481 -13.034 -2.832 12.057 31.615
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.761 6.379 2.471 0.0244 *
## strscntwt 54.303 15.408 3.524 0.0026 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.09 on 17 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.4222, Adjusted R-squared: 0.3882
## F-statistic: 12.42 on 1 and 17 DF, p-value: 0.002603
##
##
## [[3]]
##
## Call:
## lm(formula = typ3 ~ strscntwt, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.739 -18.684 -1.546 10.318 45.458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.883 9.429 2.427 0.0319 *
## strscntwt -8.900 51.381 -0.173 0.8654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.65 on 12 degrees of freedom
## Multiple R-squared: 0.002494, Adjusted R-squared: -0.08063
## F-statistic: 0.03001 on 1 and 12 DF, p-value: 0.8654